#Code built for Gabriel Fulton thesis in order to import raw text files
#Imported from text for processing, graphing, and analysis
#packages: quantmod, DT (for seperate tables)
#packages?: TTR, xts, zoo

#Set wd to whichever folder currently has the files which will be analyzed
#Surface
setwd("C:/Users/gfult/Dropbox/Salt/InProgress/Spectro Analysis/2018-09-21/PSR+3500_1596061/2018_Sep_21")
#Home PC
#setwd("C:/Users/Gabriel Fulton/Dropbox/Salt/InProgress/Spectro Analysis/2018-09-21/PSR+3500_1596061/2018_Sep_21")

#Read from .sed file based on spacing using fwf, "View(anysample)" for check
S1=read.fwf("BeetBrine_Sep21_00021.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S2=read.fwf("BeetBrine_Sep21_00022.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S3=read.fwf("BeetBrine_Sep21_00023.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S4=read.fwf("BeetBrine_Sep21_00024.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S5=read.fwf("BeetBrine_Sep21_00025.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S6=read.fwf("BeetBrine_Sep21_00026.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S7=read.fwf("BeetBrine_Sep21_00027.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S8=read.fwf("BeetBrine_Sep21_00028.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S9=read.fwf("BeetBrine_Sep21_00029.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)
S10=read.fwf("BeetBrine_Sep21_00030.sed",widths=c(-5,5,-33,10), col.names=c("Wavelength","Reflectance"),header=FALSE, skip=34)

#Graph all individual readings (Currently can't fit 10...) [layout function cool af]
#x11() makes a serpate window for graph
#x11()
#par(mfrow=c(2,5))
#plot(S1)
#plot(S2)
#plot(S3)
#plot(S4)
#plot(S5)
#plot(S6)
#plot(S7)
#plot(S8)
#plot(S9)
#plot(S10)


#Average these spectral readings with each other
purifiedwatersep21_flint_6mL=(S1+S2+S3+S4+S5+S6+S7+S8+S9+S10)/10

#Plot the cumulative average
x11()
par(mfrow=c(1,1))
plot(purifiedwatersep21_flint_6mL, main="")

#Divided by highest value such that max is now 1
reflectance=purifiedwatersep21_flint_6mL$Reflectance/max(purifiedwatersep21_flint_6mL$Reflectance)
#x11()
#par(mfrow=c(1,1))

purifiedwatersep21_flint_6mLDividedByMax=cbind.data.frame(350:2500,reflectance)
#plot()

#Correlation
#Truncating data for correlation, 350-1400
#Analysis of key points (Salt) - 1490, 1990
#this covers 350 to 1400 nm wavelength
tpurifiedwatersep21_flint_6mL=purifiedwatersep21_flint_6mL[1:1051,1:2]
cor(tsalt$Reflectance,tpurifiedwatersep21_flint_6mL$Reflectance)


#Prints to notepad to copy to excel
write.table(purifiedwatersep21_flint_6mL,"C:/Users/gfult/Desktop/Upload/purifiedwatersep21_flint_6mL.txt",sep="\t")

#Prints Normalized
write.table(purifiedwatersep21_flint_6mLDividedByMax,"C:/Users/gfult/Desktop/Upload/purifiedwatersep21_flint_6mLNORMALIZED.txt",sep="\t")



#
#
#Salt June derived POI
POIpurifiedwatersep17_5mL=rbind(salt40sep6_5mL[129,],salt40sep6_5mL[509,],salt40sep6_5mL[557,],salt40sep6_5mL[637,],salt40sep6_5mL[744,],salt40sep6_5mL[861,],salt40sep6_5mL[934,],salt40sep6_5mL[1308,])
POIpurifiedwatersep17_5mL

#ORIGINAL salt derived POI
POI2salt40sep6_5mL=rbind(salt40sep6_5mL[637,],salt40sep6_5mL[652,],salt40sep6_5mL[681,],salt40sep6_5mL[709,],salt40sep6_5mL[844,],salt40sep6_5mL[907,],salt40sep6_5mL[1108,],salt40sep6_5mL[1283,])
POI2salt40sep6_5mL

#Prints to notepad to copy to excel
write.table(salt40sep6_5mL,"C:/Users/gfult/Desktop/Upload/salt40sep6_5mL.txt",sep="\t")

#Prints Normalized
write.table(salt40sep6_5mLDividedByMax,"C:/Users/gfult/Desktop/Upload/salt40sep6_5mLNORMALIZED.txt",sep="\t")

test=rbind(salt40sep6_5mL[861,],salt40sep6_5mL[934,])
test


#Apparently embedded findPeaks is no good so adding this function:
find_peaks <- function (x, m = 3){
  shape <- diff(sign(diff(x, na.pad = FALSE)))
  pks <- sapply(which(shape < 0), FUN = function(i){
    z <- i - m + 1
    z <- ifelse(z > 0, z, 1)
    w <- i + m + 1
    w <- ifelse(w < length(x), w, length(x))
    if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])) return(i + 1) else return(numeric(0))
  })
  pks <- unlist(pks)
  pks
}
#Where m is number of points on either side of the peak
#the function can also be used to find local minima of any sequential vector x via find_peaks(-x)

#THIS AREA IS MESSED UP AND NEEDS WORK
#Peak and valley analysis (NOT SURE OF WHAT m IS APPROPRIATE, originally 3)
#peaks=find_peaks(Stotal$Reflectance, m=30)
#valleys=find_peaks(-Stotal$Reflectance, m=30)

#peaks=findPeaks(Stotal$Reflectance,thresh=0.05)
#valleys=findValleys(Stotal$Reflectance,thresh=0.05)

#Finds corresponding reflectance values for found wavelength peaks and valleys
#ppeaks=Stotal[Stotal$Wavelength%in%peaks,]
#pvalleys=Stotal[Stotal$Wavelength%in%valleys,]


#Highlight of said peak and valleys with corresponding wavelengths/severity
#points(ppeaks,col="red",pch=19)
#points(pvalleys,col="blue",pch=19)

#Creation of tables of corresponding data:
#Peaks:
#View(ppeaks, "Peaks")
#Valleys:
#View(pvalleys, "Valleys")

#For Table writing:
#write.xlsx(x, file, sheetName="Sheet1")

